Cis-Regulatory Modules Drive Dynamic Patterns of a Multicellular System 
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How intracellular and extracellular signals are integrated by transcription factors is essential for 
understanding complex cellular patterns at the population level. In this Letter, by using a synthetic 
genetic oscillator coupled to a quorum-sensing apparatus, we propose an experimentally feasible 
cis-regulatory module (CRM) which performs four possible logic operations (ANDN, ORN, NOR 
and NAND) of input signals. We show both numerically and theoretically that these different CRMs 
drive fundamentally different dynamic patterns, such as synchronization, clustering and splay state. 
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Biological organisms possess an enormous repertoire of 
genetic responses to ever- varying combinations of cellular 
and environmental signals Such a repertoire is typ- 

ically encoded in complex regulatory networks, and af- 
fects patterning, differentiation and growth. At the heart 
of these networks are cis-regulatory modules (CRMs), 
which contain a cluster of binding sites for transcrip- 
tion factors (TFs) and determine the place and timing 
of gene action within the network. Both deciphering the 
codes and elucidating the functions of CRMs involved in 
various developmental processes are a major challenge in 
biology. 

It has been shown that CRMs can perform an elaborate 
computation at the individual gene level: the transcrip- 
tion rate of a gene depends on the active concentration of 
each of inputs 3, 4, 5, 6, 7]. On the other hand, cells live 
in a complex environment and can sense many different 
signals, in particular those from neighboring cells. There- 
fore, at the multicell level CRMs need to integrate intra- 
cellular and extracellular signals so as to coordinate gene 
expression. Given that cells are frequently subject to 
chemical signals from neighboring cells, it is worth study- 
ing the effect of chemical communication on the dynamic 
patterns of multicellular systems. Modeling studies, for 
example, have shown that a population of rcpressilators 
coupled to quorum sensing can work as a macroscopic ge- 
netic clocks In that study, two input signals (i.e., two 
TFs) regulate a target gene independently. TFs, how- 
ever, are often integrated in a combinatorial logic man- 
ner, and moreover such a combination may take different 
forms [3, 4, 5, 6, 7]. From views of evolutionism, CRMs 
are changeable, e.g., cis-regulatory mutations 0. Such 
a mutation constitutes an important part of the genetic 
basis for adaptation. A naturally arising question is how 
the changes of CRMs affect cellular patterns of popu- 
lations of genetic oscillators. We address this question 
by designing a multicellular network with a CRM con- 
sisting of repressilators coupled to quorum sensing 
[^, IllLfl^ in Escherichia coli. In contrast to the previous 



studies [i, HI that numerically showed that coupled ge- 
netic oscillators can demonstrate synchronous behaviors, 
we both numerically and theoretically show that differ- 
ent signal integration (ANDN, ORN, NOR, NAND type 
of responses) leads to fundamentally different properties, 
such as synchronization, clustering, and splay state. Our 
results indicate that the CRM has a significant influence 
on the mode of cellular patterns. 

A multicellular network under investigation is 
schematically shown in Fig. 1(a). In such a network, 
the signaling molecule (S) carries out the information 
exchange between cells and regulates the expression of a 
target gene through a CRM. The S and the TF (Y) first 
bind to specific DNA sequences of the CRM, and then 
co-regulate the expression of the gene in a combinatorial 
scheme. In theory, this type of CRM can perform eight 
different cis-regulatory input functions (CRIFs) ^^], but 
limited by the cyclic repression structure of repressilator, 
we have only four types of CRIFs: ANDN, ORN, NOR 
and NAND (see Ref. for exact explanations). Figure 
1 (b) gives the detailed regulation scheme of every CRM. 

Based on the biochemical reactions given in Table 1 
and defining the rescaled concentrations as our dynami- 
cal variables, the dimensionless equations of intracellular 
dynamics are described as 

= F(X„5,), (1) 



dt 
~dt 



(2) 



where subscript 



1,2,- ■• ,iV), and 
Ui and Zi stand- 



i represents cell i (= 
Xi, Yi, Zi)^ with X. 
ing for three mRNA concentrations, and X^, and Z.^ 
for three protein concentrations. Si represents the con- 
centration of the signaling molecule inside the ith cell 
whereas Se does the concentration of the signal in the 
extracellular environment. Because of the fast diffu- 
sion of the extracellular signal compared to the repres- 
silator period, Se can be assumed to be in the quasi- 
steady state, leading to Se — where the pa- 



2 




RNAPol 



RNAPol RNAPol.-'"" 








P, o„ 

RNAPol 






V 1 


n 1 


P o„, 0„, 

RNAPol 

(SJj (jjs: 




1 m 


n 





(b) 

ANDN 

ORN 
NOR 

NAND 



FIG. 1: (color), (a) The schematic diagram of a multicellular system with a cis-regulatory module (CRM). Three transcriptional 
repressors (X, Y, Z) inhibit one another in a cyclic way. The gene /us/ from the LuxI/LuxR module first synthesizes a small 
molecule S. Both S and LuxR then form a hetero-tetramer complex. The dimer of Y and the complex co-regulate the target 
gene, thus carrying out the function of a CRM (symboUed by the empty box). The right bidirectional arrow indicates that S 
can freely diffuse through the cellular membrane, (b) Four cis-regulatory constructs for implementations of four different logic 
functions. From top to bottom is ANDN, ORN, NOR and NAND, respectively. In ANDN, the activator S is unable to act if 
the repressor Y is bound to the promoter; In ORN, the CRM is constructed using a weak promoter and a strong promoter, 
where the activator S is able to act if the repressor Y is bound to the promoter; In NOR, two repressors S and Y produce the 
full repression cooperatively; In NAND, the promoter is regulated exclusively by two repressors. In (a) and (b), P and Or 
denote the promoter (dark blue box for the strong promoter and light blue box for the weak promoter) and the operator site 
(jacinth box) respectively, and the RNA Pol represents RNA polymerase. We use offset and overlapping boxes to indicate the 
mutual repression and the dashed lines to indicate the cooperative interaction. 



rameter Q depends on the cell density in a nonlinear 
way i. F = (Fi,F2,...,F6)T with Fi = jf^ - x, 

F2 = -y, Fi ^ f5{x - X), Fs = f3{y - Y), 

Fe = I3{z ~ Z), E ^jX ~ 6S, and 



Fa = CRIF 



(3) 



where we omit subscript i for convenience, and the core 
function CRIF corresponding to ANDN, ORN, NOR and 
NAND respectively is listed in the first part of Table 
1. The detailed derivation of CRIFs and F is put in 
Ref. Throughout this Letter, all parameters except 
for Q are set as a = 204, /? = 1, 7 = 0.01, 5=1, 
n = 2, ?7 = 2, /i = 51, v = 204, A = 1, which come 
from experimentally-reasonable settings [13]. Since the 
numerical results do not depend qualitatively on the cell 
number, we set TV = 120. 

We are interested in the influence of four possible 
CRIFs on cellular patterns. The results shown in the 
insets of Fig. 2 indicate that these different CRMs drive 
fundamentally different dynamic patterns. Specifically, 
in the case of ANDN, for arbitrarily chosen initial con- 
ditions we observe complete synchronization (1-cluster) 
only, similar to that shown in Refs. [1, [HI, This 
pattern indicates that a specific CRM would combine in- 
tracellular and intercellular signals to coordinate the gene 
expressions in a uniform way at the population level. In- 
terestingly in the case of ORN, we find that different ini- 
tial conditions lead to three kinds of dynamic patterns: 
1-cluster, 2-cluster and 3-cluster 15|. Similar phenomena 
were also found in a chemical system [l^, [l3| . In the case 



TABLE I: Biochemical reactions and cis-regulatory input 
functions (CRIFs). See Ref. ^ for the derivation of CRIFs, 
experimental values of parameters (including /i, 1/ and A that 
depend on reaction rates), and meanings of all used symbols. 
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due to different promotors. 
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Kuramoto phase reduction method 2l| gives 



FIG. 2: (color). Different cis-regulatory modules drive dif- 
ferent cellular patterns. Insets display instantaneous distri- 
butions of phases of the oscillators for a fixed Q — 0.5: (a) 
1-cluster state (complete synchronization) for ANDN; (b) 1- 
, 2- and 3-cluster states for ORN; (c) splay state for NOR; 
(d) 3-, 4- and 5-cluster states for NAND (the corresponding 
time courses of all clustering figures are put in the support- 
ing material Ref. TS]). The function G{A(j}) determines the 
coupling mode: attractive coupling for ANDN and ORN (due 
to G'(0) > 0) and repulsive coupling for NOR and NAND 
(due to G'{0) < 0). Here, the different clusterings arise from 
different initial phases. 



of NOR, however, neither synchronization nor clustering 
is observed, but an interesting phenomenon that all cells 
are staggered equally in time, i.e., so-called splay state, is 
found for the first time in a cell population although the 
similar phenomenon was also detected experimentally in 
a multimode laser system Finally, in the case of 

NAND, we also observe three types of clusterings at the 
scattered initial states: 3-, 4- and 5-clusters 3, 23) ■ The 
complete synchronization, however, never occurs in this 
case. 

To understand and interpret the above interesting pat- 
terns, we have performed an analytical study of the sys- 
tem in the phase model description [iH , which holds in a 
weak coupling case. In this description, we first rewrite 
Eq. (2) as the following symmetric form of coupling 



= E{X,,S,) - 77(1 - Q)S. + - ^ rjQiS, - S.) . (4) 



Then, for convenience the system consisting of both Eq. 
(1) and the equation 



at 



(5) 



dt 



1 ^ 



(6) 



is called as auxiliary system, which is assumed to gen- 
erate a sustained oscillation. For a weak coupling, the 



where and uji stand for the phase and frequency of the 
auxiliary system, respectively. iJij(A0) represents the 
interaction function with respect to the phase difference 
A(j) = 4>j — (pi between two cells, 

H,, (0j _ 0,) = _ / z(e) ■ p{(f>, - 0, + e)de (7) 

which can be calculated numerically [22], where Z{d), 
a phase response function characterizing the phase ad- 
vance per unit perturbation, is a 27r-period function, and 
p = (0,0,0,0,0,0,r7Q(S'j - S',))'^. Below we will omit 
subscripts i and j for convenience. From H{A(j)), we 
introduce a function: G(A0) = -ff(A0) - i7(-A0), to 
determine the mode of coupling. If G{A(j)) exhibits a 
positive slope at Ap = 0, i.e, G'(0) > 0, the couphng 
is phase-attractive; If G'(0) < 0, the coupling is phase- 
repulsive. Therefore, Fig. 2 implies that the CRMs in 
the cases of ANDN and ORN correspond to the phase- 
attractive coupling whereas those in the cases of NOR 
and NAND correspond to the phase-repulsive coupling. 
Such an approach based on the sign of G'(0) that depends 
generally on the intrinsic dynamics of the uncoupled os- 
cillator and on the interaction between the oscillators is 
more effective than that of directly observing the network 
topology in determining the mode of weak coupling [l^ , 
especially in the case of complex network architectures. 

One cannot, however, obtain knowledge about cluster- 
ing from the sign of G'fO). Since we are interested mainly 
in balanced clusters [l5|, we next employ Okuda's ap- 
proach to determine the stability of such clusters [1^ . In 
that method, we need to calculate two kinds of eigen- 
values: one is associated with intra-cluster fluctuations 
and the other with inter-cluster fluctuations, which are 
denoted by Xp and Xq (see the caption of Fig. 3) re- 
spectively, where M < p < N — 1 and < q < M — 1 
with M being the number of clusters presumptively. For 
convenience, denote by A^^^ and A^^) the iV - M same 
eigenvalues Ap and the maximum of the the real parts of 
(M — 1) non-zero eigenvalues Xq, respectively. Then, the 
stability of clusterings can be determined by the signs of 
A'-^-' and A'^' . Specifically, the clustering is stable if both 
A'-^-' and A'-^-' are negative, and unstable if A^^^ is positive. 
In addition, if A^^^ is positive and A*^^-' is negative, and 
further if M ~ N, the M-cluster (i.e., the splay state) 
are also stable. The dependence of A*^^^ and A*^^^ on the 
balanced cluster number M in the cases of four CRIFs is 
shown in Fig. 3, which further verifies the dynamic pat- 
terns shown in Fig. 2. The insets of Fig. 3 show that the 
parameter Q has the significant influence on the stability 
of clusterings (even including 1-cluster in Fig. 2(a) and 
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terns and natural developmental processes. 
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FIG. 3: (color). Eigenvalues associated with intra-cluster 
fluctuations (A'^-*: blue circle) and the maximal real part of 
non-zero eigenvalues associated with inter-cluster fluctuations 
(A'^'': red square) as a function of the number of balanced 
clusters (M) for a fixed Q = 0.5 in the case of: (a) ANDN; 
(b) ORN; (c) NOR and (d) NAND. Insets: the dependence 
relation of A'^-* and A'^' on the parameter (Q) for a particular 
clustering as indicated. 



the splay state in Fig. 2(c)) for a particular balanced 
cluster state, according to the above analysis. 

In addition, in order to verify that the above results 
are of generality, we also investigated the case of genetic 
relaxation oscillators by using a detailed example studied 
in Ref. , and found that different CRMs also drive fun- 
damentally different dynamic patters, but different types 
of CRIFs would lead to different cellular patterns from 
those in the case of repressilator (due to the paper length, 
the detailed results are displayed in [13]). 

In summary, using models of synthetic genetic oscil- 
lators coupled to quorum sensing, we have shown that 
different CRMs drive fundamentally different cellular 
patterns, such as synchronization, clustering, and splay 
state. Our results imply the following two points: (1) 
Multicellular organisms possibly evolve into some func- 
tional CRMs for particular goals by performing an elab- 
orate computation for input TFs; (2) Genetic network 
architecture found in synchronous circadian clocks [1, [^l 
might be constrained since the complete synchroniza- 
tion independent of initial conditions takes place only 
in the case of ANDN type of responses. In particular, 
our results do suggest possible candidate circuits for syn- 
chronous circadian clocks, while excluding others. We 
expect that our theoretical findings will stimulate further 
investigations under a more realistic condition involving 
stochasticity [25', '26'] and spatial heterogeneousness 27 1 , 
which would help us to understand differentiation pat- 
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